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Abstract 

In galactic nuclei with sufficiently short relaxation times, binary supermassive black holes can evolve beyond 
their stalling radii via continued interaction with stars. We study this "collisional" evolutionary regime using 
both fully self-consistent A^-body integrations and approximate Fokker-Planck models. The A^-body integra- 
tions employ particle numbers up to 0.26 x 10*' and a direct-summation potential solver; close interactions 
involving the binary are treated using a new implementation of the Mikkola-Aarseth chain regularization al- 
gorithm. Even at these large values of A^, two-body scattering occurs at high enough rates in the A^-body 
simulations that the binary is never fully in the diffusively-repopulated (i.e. large-A^) loss cone regime, which 
precludes a simple scaling of the results to real galaxies. The Fokker-Planck model is used to bridge this gap; 
it includes, for the first time in this context, binary-induced changes in the stellar density and potential. The 
Fokker-Planck model is shown to accurately reproduce the results of the A^-body integrations, and is then ex- 
tended to the much larger A^ regime of real galaxies. Analytic expressions are derived that accurately reproduce 
the time dependence of the binary semi-major axis as predicted by the Fokker-Planck model. Gravitational 
radiation begins to dominate the binary's evolution after a time that is always comparable to, or less than, the 
relaxation time measured at the binary's gravitational influence radius; the observed correlation of nuclear re- 
laxation time with velocity dispersion implies that coalescence in < 10 Gyr will occur in nuclei with a 80 
km s^', i.e. with binary black hole mass <i2x lO^Mp^. The coalescence time depends only weakly on binary 
mass ratio. Formation of a core, or "mass deficit," is shown to result from a competition between ejection 
of stars by the binary and re-supply of depleted orbits via two-body scattering. Mass deficits as large as ~ 4 
times the binary mass are produced before the gravitational radiation regime is reached; however, after the 
two black holes coalesce, a Bahcall-Wolf cusp appears around the single hole in approximately one relaxation 
time, resulting in a nuclear density profile consisting of a flat core with an inner, compact cluster, similar to 
what is observed at the centers of low-luminosity elliptical galaxies. We critically evaluate recent claims that 
binary-star interactions can induce rapid coalescence of binary supermassive black holes even in the absence 
of loss cone refilling. 
Subject headings: 



I. INTRODUCTION 

This paper is the third in a series investigating the evolution 
of binary supermassive black holes at the centers of galaxies. 
A massive binary hardens via exchange of energy and angular 
momentum with passing stars, but this process is self-limiting, 
since the interacting stars are ejected from the nucleus with 
velocities of order the relative velocity of the two black holes. 
Continued hardening of the bi nary requires a repopulation of 
the depleted orbits. Paper I CMilosavljevic & Mertt3l2003h 
discussed various mechanisms by which this can occur, in- 
cluding collisional loss-cone repopulation, secondary sling- 
shot, chaotic stellar orbits, and Brownian motion of the bi- 
nary. These different mechanisms typically obey different 
scalings of the binary hardening rate with the number A^ of 
stars and with time; in the large-A^ limit and in a spherical or 
axisymmetric potential, the hardening rate (defined as the rate 



of change of the binary's energy) is predicted to scale roughly 
as A^^ ' , i.e. inversely with the relaxation time, and hence to be 
very small for values of A' characteristic of massive elliptical 
galaxies ( Valtonen 1996; Yu 2002). 

In Paper II ( Berczik et al.l l2005l) . a direct-summation A^- 
body code, combined with a parallel GRAPE cluster, was 
used to carry out integrations of binary evolution in galaxy 
models with large, low-density cores. Because of their low 
central density, the relaxation time at the center of these mod- 
els was relatively long (compared with orbital periods), and 
collisional loss cone refilling was shown to occur at a lower 
rate than the loss of stars to the binary, i.e. the binary's loss 
cone remained nearly empty. This is the same ("diffusive") 
regirne believed to characterize binary evolution in real galax- 
ies (iMilosavlievic & Merritli 120011) . The A^-body hardening 
rates were compared with the predictions of simple loss-cone 
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theory and found to be in reasonable agreement. 

The Plummer models used in Paper II were not good rep- 
resentations of real galaxies. In this paper, we present a new 
set of simulations based on galaxy models that more closely 
approximate real galaxies, with power-law central density 
cusps. In order to deal efficiently with interactions involv- 
ing the binary, we incorp orate the Mikkola-Aarseth chain- 
regularization algorithm jMikkola & AarsethI 1 19901 Il993h 
into our A^-body code, including both the effects of nearby 
stars as perturbers of the chain, and the effects of the chain 
on the surrounding stars. The resulting A^-body algorithm 
is coupled with a GRAPE-6 special-purpose computer and 
used to carry out extended integrations of binaries with var- 
ious values of A^, up to the limit « 0.26 x 10*' set by the 
grape's memory. In order to more accurately characterize 
the A^-dependence of the evolution, multiple integrations are 
carried out starting from different random realizations of the 
same initial conditions and averaged. 

Even at the large values of A^ allowed by the GRAPE-6, 
two-body (star-star) scattering occurs at a high enough rate in 
the A^-body simulations that the binary is never fully in the 
empty-loss-cone regime. This fact precludes a simple scaling 
of the A^-body results to real galaxies. We therefore develop 
a Fokker-Planck model that can be applied to nuclei with any 
value of A^, i.e. any value of Mi2/»J*, where Mn = Mi +M2 
and are the mass of the binary and of a single star, respec- 
tively. Our Fokker-Planck model is unique in that it allows 
for the joint evolution of the binary and of the stellar nucleus; 
it can t herefore reproduce the cre ation of a core, or "mass 
deficit" (Milos avlievic et aT]|2002h . as the binary ejects stars. 
The Fokker-Planck model is first tested by comparison with 
the A^-body results, and is then applied to the much larger-A^ 
regime of real galaxies. In this way we are able to make the 
first detailed predictions about the joint evolution of massive 
binaries and stars at the centers of galaxies. 

The time scale that limits binary evolution in our models 
is the relaxation time, defined as the time for (mostly distant) 
gravitational encounters between stars to establish a locally 
Maxwellian velocity distribution. Assuming a homogenous 
and isotropic distribution of equal-mass stars, the relaxation 
time is approximately 
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(ISpitzeill98 7). Here, Oioo is the Id stellar velocity dispersion 
in units of 100 km s^\ ps is the stellar mass density in units 
of IO^Mq pc^^, OT* = mi,/MQ, and InAis = In A/ 15, whe re 
InA is the Coulomb logarithm and A w 0.4A^ (ISpitzedll987h . 

Figure [U shows estimates of T^, measured at the supermas- 
sive black hole's influence radius r/,, for the ACSA^irgo sam- 
ple of early-type galaxies (ICdte et al.ll2004h . The influence 
radius was defined in the usual way via 



M*(r,0-2M. 



(2) 



and the black hole mass was inferred from the measured value 
of o via the the M, — o relation, 

M. w 5.72 X IO^MqcH^ (3) 

(iFerrarese & FordI l2005h . A stellar mass of IMq was as- 
sumed. 

Figure [T]reveals a tight correlation between ri (r/,) and a. A 
least-squares fit to the points (shown as the dashed line in the 
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Fig. 1 . — Relaxation times, measured at the supermassive black hole's in- 
fluence radius, in the ACSA'irgo sample of galaxies (Cote et al. 2004), versus 
the central stellar velocity dispersion. Filled symbols are galaxies in which 
the black hole's influence radius is resolved; star is the Milky Way. 



figure) gives 



r,(r,) « 1.16 X 10^1 yraloo' 
«8.0x 10''yrM^-f 



(4a) 
(4b) 



where M, g = M,/ IO^Mq. The results presented in this paper 
are only relevant to galaxies in which the nuclear relaxation 
time is not much longer than galaxy lifetimes; according to 
Figure [1] this is the case for galaxies with a ^ 80 km s^'. 
This is roughly the velocity dispersion near the center of the 
Milky Way; hence, the sort of evolution that is modelled here 
is most relevant to spheroids that are not much brighter than 
the Milky Way bulge. 

The A^-body techniques are described in §2 and §3 and the 
results of the A^-body integrations are presented in §4. In §5 
the Fokker-Planck model is described and compared with the 
A^-body results. Predictions of the Fokker-Planck model in the 
large-A^ regime corresponding to real galaxies are presented in 
§6. §7 and §8 discuss the implications for evolution of binary 
supermassive black holes in real galaxies, and §9 presents a 
critical comparison with other proposed models of binary evo- 
lution. §10 sums up. 

2. W-BODY TECHNIQUES 

Our A^ -body algorithrn was an adaptation of the NBODYl 
code of I AarsethI (Il999h to the GRAPE-6 special purpose 
hardware. The code uses a fourth-order Hermite int egration 
scheme with individual, adaptive, block time steps dAarsethl 
l2003h . For the majority of the particles, the forces and force 
derivatives were calculated via a direct-summation scheme 
using the GRAPE-6. More details of the particle advance- 
ment scheme can be found in Paper II. As discussed there, the 
code contains two parameters that affect the speed and accu- 
racy of the calculation, the particle softening length e and the 
time-step accuracy parameter r|. 

Close encounters between the massive particles ("black 
holes"), or between black holes and stars, require pro- 
hibitively small time steps in such a scheme. To avoid this 
situation, we adopte d a chain re gularization algorithm for the 
critical interactions tMikkola & Aarseth 1990„ 1993) , as fol- 
lows. Let r,, / = 1,...,N be the position vectors of the par- 



Binary Black Holes 



3 



tides. We first identify the subset of n particles to be in- 
cluded in the chain; the precise criterion for inclusion is pre- 
sented below, but in the late stages of evolution, the chain 
always included the two black holes as its lowest members. 
We then search for the particle that is closest to either end of 
the chain and add it; this operation is repeated recursively un- 
til all n particles are included. Define the separation vectors 
R, = r,+ i — r, where r,+i and r,- are the coordinates of the two 
particles making up the /th link of the chain. The canonical 
momenta W,- corresponding to the coordinates R,- are given in 
terms of the old momenta via the generating function 



n-l 



IW,.(: 



(5) 



Next, we apply KS regularization (iKustaanheimo & Stiefell 
Il965h to the chain vectors, regularizing only the interactions 
between neighboring particles in the chain. Let Q, and P,- be 
the KS transformed R, and W, coordinates. After applying 
the time transormation 5f = ghs, g ^ l/L, where L is the La- 
grangian of the system {L = T — U, where T is the kinetic and 
U is the potential energy of the system). We obtain the reg- 
ularized Hamiltonian F = §(//(Q,,P,) —Eq), where Eq is the 
total energy of the system. The equations of motion are then 



(6) 



where primes denote differentiation with respect to the time 
coordinate s. Because of the use of regularized coordinates, 
these equations do not suffer from singularities, as long as 
care is taken in the construction of the chain. 

Since it is impractical to include all particles in the chain, 
we must consider the effects of external forces on the chain 
members. Let Fj be the perturbing acceleration acting on the 
y'th body of mass nij. The perturbed system can be written in 
Hamiltonian form by simply adding the perturbing potential: 



(7) 



Only one chain was defined at any given time. At the start 
of the A^-body integrations, there was no regularization, and 
all particles were advanced using the variable-time-step Her- 
mite scheme. The first condition that needed to be met before 
"turning on" the chain was that one of the particles (including 
possibly a black hole) achieved a time step shorter than ?£■/„„,„ 
and passed a distance from one of the black holes smaller than 
rdimin- If this Condition was satisfied, it was then checked 
whether the encounter resulted in a deflection angle greater 
than 25 = 7t/2, where 



cos 5 = 



1 



-1/2 



(8) 



here R is the impact parameter, Vq is the pre-encounter relative 
velocity, and m\ and m2 are the masses of the two particles. 
This condition is equivalent to 



mi +11X2 >RVq. 



(9) 



Each star closer to the black hole than rchmin was then added 
to the chain, and the two black holes were always included. 
The values of tchmin and rchmin were determined by carrying 
out test runs; we adopted fe/imm ~ 10^^ — 10^^ and rchmin ~ 
10^^ — 10^^ in standard A^-body units. 
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Fig. 2. — The average time step, as defined in the text, during two in- 
tegrations of a binary l^lack hole at the center of a Dehen-model galaxy. 
N = 20,000, and the softening length and time-step parameters of the N- 
body code were E = 10^*, T| = 0.01. In the absence of the chain, the average 
time step drops to very low values once the binary begins to harden. 
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Fig. 3. — Relative energy error over 100 time units of of a set of integrations 
like those in Fig. |2] for various values of the softening length e, and with the 
chain. 



The chain's center of mass was a pseudoparticle as seen by 
the A^-body code and was advanced by the Hermite scheme in 
the same way as an ordinary particle. However, when inte- 
grating the trajectories of stars near to the chain, it is essential 
to resolve the inner structure of the chain. Thus for stars inside 
a critical rc-iti radius around the chain, the forces from the in- 
dividual chain members were taken into account. The value 
of rcriti was set by the size of the chain to be rcriti — XRch 
with Rch the spatial size of the chain and X = 100. In addi- 
tion, the equations of motion of the chain particles must in- 
clude the forces exerted by a set of external perturber stars. 
Whether or not a given star was listed as a perturber was deter- 
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mined by a tidal criterion: r < Rcriti = {m/mchainY^^lJil^ Rch 
where nichain represents the mass of the chain, m is the mass 
of the star, and Y„„>, was chosen to be 10^^; thus rcrui ~ 

The membership of the chain changed under the evolution 
of the system. Stars were captured into the chain if their or- 
bits approached the binary closer than 7?^/,. Stars were emitted 
from the chain if they got further from both of the black holes 
than 1 -SRch- The difference between the emission and absorp- 
tion distances was chosen to avoid a too-frequent variation of 
the chain membership. When the last particle left the chain, 
the chain was eliminated and the integration turned back to 
the Hermite scheme, until a new chain was created. 

In what follows, we refer to the A^-body code without chain 
as NBl, and the code including chain as CHNBl. We carried 
out a number of tests to see how the performance and accu- 
racy of the NB 1 code were affected by inclusion of the chain. 
Typically, one integration step of the chain required about 
five times as much cpu time as a single call to the GRAPE- 
6, due to the complex nature of the chain and the generally 
large number of perturber particles. Thus our code is quicker 
than a basic Hermite scheme code (NBl) only if the small- 
est time steps are about an order of magnitude shorter than 
the next smallest time steps, and if in addition those particles 
would be assigned to the chain. It is easy to show that in the 
case of a galaxy including a central, massive binary system 
this condition is usually fulfilled. The Hermite time step of 
the binary is considerably smaller than the time steps of the 
stars, due to their close orbit and fast evolution. Of course, 
the performance of both codes depends on the two parameters 
r| (time step parameter) and £ (particle softening length) that 
determine the accuracy of the star-star interactions. In what 
follows, we fixed r| = 0.01 based on the results of the tests in 
Paper II. In CHNBl, e was always set to zero. 

Figures |2l|4] show the results of our performance tests. Fig- 
ure [2] plots the average time step as a function of time in 
both codes, for integrations of a binary black hole in a galaxy 
model following Dehnen's (1993) density law: 
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with y = 1 .2 and N = 20, 000 particles. The two black holes 
had equal masses. Mi = Mi = 0.01M|^,a/, and were placed in- 
tially on a circular orbit with separation 0. lOa. We defined the 
average time step as t/Nnmeaeps (t), where Ntimesteps (0 was the 
total number of integration time steps until time f, including 
only the time steps of particles outside the chain. It can be 
seen that in the early stages of the evolution, the time steps are 
about the same in both cases. However as the binary becomes 
harder, the NBl time steps become smaller and smaller, in 
order to achieve the necessary precision in the integration of 
the binary. In the code with the chain, the average time step 
hardly changes after the binary begins to harden. The binary 
is integrated by the regularized equations, hence the step size 
of the NBl integration remains relatively large. The average 
step size of the CHNBl code was about 2.5 x 10^^ (in units 
where G = Mf,ai = a = 1), while in the NBl integration it was 
3 X 10^^. The resulting net speed-up with the chain is more 
than a factor of two. 

Figure [3] shows the energy error in a set of integrations of 
the same model but with various different softening lengths 
e, compared with the energy error in an integration with the 
chain (and e = 0). It can be seen that the energy conservation 
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Fig. 4. — Results of a set of test integrations with and witliout the chain. 
Initial conditions consisted of a binary of mass M\ = Mi = 0.005 and sepa- 
ration 0. 1 , in a Dehnen-model galaxy with y = 1 .2 and N = 20, 000. 



of the CHNBl code is about as good as the best results from 
the NBl integrations. However, the latter occur when the soft- 
ening length is very large, much too large for accurate integra- 
tion of the binary. This is shown in FigureH] It is evident from 
that figure that with larger softening lengths, e ^ 10^^, the in- 
tegration of the binary is not very accurate. However, even 
with very small e, the evolution of the distance between the 
black holes includes "spiky" behaviors, due apparently to the 
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Fig. 5. 
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- Evolution of 1/a, the inverse binary semi-major axis (left column), and e, the binary eccentricity (right column), in the full set of A'-body integrations. 
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very sensitive nature of the eccentricity evolution with respect 
to the precision of the integration (Figure|4]3,c). 

These results suggest that chain regularization is an accu- 
rate and efficient way to integrate binary black holes at the 
centers of galaxies. It can keep track of the evolution of the 
binary with high precision, and the calculation time is sub- 
stantially faster than a plain Hermite integration when the lat- 
ter is used with a reasonable (i.e. sufficiently small) softening 
parameter. 

3. INITIAL CONDITIONS 

All of our integrations adopted Dehnen's model, equa- 
tion (fTOb . for the initial galaxy, with y = 0.5. To this model 
were added two particles, the "black holes," with masses 
Ml = M2 = 0.005 in units of the galaxy mass. (Henceforth 
we write M12 = Mi +M2.) The black holes were placed 
symmetrically about the center of the galaxy at x = ±0.1. 
The initial velocities of the black holes were chosen to be 
Vy = ±0.16 yielding nearly circular initial orbits. These initial 
conditions are similar to those adopted in some earlier studies 
dOuinlan & Hernquist ( 1997); Nakano & Makino ( 1999)) al- 
though they are probably less realistic than initial conditions 
that plac e one of the two massive particles exactly at the cen- 
ter (e.g. iMerritt & Szelll (120061) ). Henceforth we adopt units 
such that the gravitational constant G, the total mass in stars 
Mga/, and the Dehnen scale length a are equal to one. In these 

units, the crossing time {GMg,,! / a^y^/^ is also equal to one. 

A standard expression for r/,, the radius of influence of a 
single black hole at the center of a galaxy, is 

M*(r,,)=2M. (11) 

with M, the black hole mass and Mj,(r) the mass in stars 
within a sphere of radius r . The semi-major axis length of 
a "hard" binary is som etimes defined in terms of r/, as (e.g. 
IMerritt & WangI 120051) ) 

a r/, 



TABLE 1 

Parameters of the A'-body integrations 



ah 



(1+a) 



(12) 



with a = M2/M1 < 1 the binary mass ratio, and we adopt that 
definition here. SettingM. =Mi2 = 0.005 +0.005 = 0.01 and 
a = 1 , the values of r/, and a/, for our A^-body models are 

r/, =0.264, fl/, =0.0165. (13) 

These expressions ignore the changes that the two black holes 
induce in the mass distribution of the galaxy when forming a 
hard binary, but are useful as points of reference. 

Based on the results of Papers 1 and 11, once the binary has 
interacted with and ejected most of the stars on intersecting 
orbits, its subsequent evolution is dependent on the contin- 
ued scattering of stars into its sphere of influence; since the 
scattering time scale increases with A^, the binary's decay rate 
should decrease as increases. In order to better characterize 
this A^-dependence, the initial conditions were realized using 
six different values of N, N = (8192, 16384, 32768, 65536, 
131072, 262144), otN = 2P, p= (13, 14, 15, 16, 17, 18). (In 
what follows, we refer to these different A^-values via the 
shorthand 8K, 16K, 262K). The lai-gest of these A^ val- 
ues is close to the maximum number of particles that can be 
handled in the GRAPE-6 memory. In order to decrease the 
"noise" associated with the evolution for small A^, we carried 
out n,„, multiple integrations at each A^, in which the initial 
stellar postions and velocities were calculated using different 
seeds for the random number generator All of these integra- 
tions were continued until a time T,„ax = 350; when scaled to 



Name 


N 


"»» 


8K 


8192 


18 


16K 


16384 


8 


32K 


32768 


6 


65K 


65536 


4 


131K 


131072 


2 


262K 


262140 


1 



a typical luminous elliptical galaxy with crossing time 10^ 
yr, this corresponds to ^ lO"' yr. Table 1 gives the parameters 
of the A^-body integrations. 

4. W-BODY RESULTS 

Initially the two black holes move on nearly independent 
orbits about the center of the galaxy. The orbits decay, and 
at f w 10 the black holes form a hard binary. After this, the 
semimajor axis a of the binary shrinks as the two black holes 
interact with stars and eject them from the nucleus via the 
gravitational slingshot. Figure |5] shows the evolution of l/a 
and e, the orbital eccentricity, in the full set of integrations for 
t > 20. The scatter in the values of l/a and e at a given time 
is considerable in the integrations with smallest A^. Neverthe- 
less a clear trend is apparent: both l/a and e evolve less, on 
average, as A^ is increased. 

4. 1 . Binary Hardening 

In order to clarify the A^-dependence of the evolution, we 
computed averages over the n,„, independent integrations of 
fl^'(f) and e(f). Figure |6] shows the mean evolution of l/a 
for the six different A^ values. The early evolution (Fig. |6p), 
until f 10, is essentially A^-independent. In this regime, the 
hardening of the binary is driven by dynamical friction against 
the stars, and the rate of binding energy increase is a function 
only of the stellar density, which is the same for each of the 
A^-body models. 

At f « 10, the binary hardening rate begins to show a clear 
A^-dependence, in the sense of more gradual hardening for 
larger A^. In Merritt (2006), the separation at which this oc- 
curs was defined as the "stalling radius," since in the limit of 
large A^ the binary would stop evolving at this point. Based on 
Figure|6l aj^l,, w a,;' w 60 and t„„ii w 10. 

At f i> tstaii the A^-dependence of the evolution is striking 
(Fig.lJb). As in Paper II, we define the instantaneous harden- 
ing rate as 

Figure I2] shows {s){t) computed by fitting smoothing splines 
to the averaged a^^{t) curves of Figure [6)3. Mean harden- 
ing rates are roughly constant with time for each A^. The A^- 
dependence of the hardening rate is shown in Figure [8] Here, 
(5) was computed by fitting a straight line to (fl^')(f) in an 
interval Af = 50 centered on (a^')(r) = 750; in this way, the 
different hardening rates are being compared at similar val- 
ues of the binary semi-major axis, chosen to be roughly the 
minimum value reached in the integration with largest A^. The 
dependence of (s) on A^ is approximately a power law. 



log 



10-' 



2.27- 0.357 logioA^. 



(15) 



The ~ A^ '^■^ dependence is considerably flatter than the ~ 
A'^' dependence expected in a diffusively -refilled loss cone 
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Fig. 6. — Short-term (a) and long-term (b) evolution of the mean value of 1/a in the W-body integrations. Horizontal line in panel (a) indicates approximately 
where the transition occurs between W-independent and W-dependent evolution; this is also roughly the "stalling radius" defined in Merritt (2006), and the "hard 
binary" separation defined in Yu (2002). 



dMilosavlievic & Me rritt 2003). This fact precludes any sim- 
ple extrapolation of the data in Figure [8] to the much larger 
regime of real galaxies. 

We can compare these hardening rates with the predictions 
of scattering experiments in a fixed, infinite, homogeneous 
background: 



d 
dt 



H 



Gp 



(16) 



with p and G the mass density and Id velocity dispersion of 
the stars, and H a dimensionless rate coefficient that depends 
on the binary separation, mass ratio and eccentricity. For a 
hard, equal-mass, circular- o rbit binary, H » 16 {Hills 19831 
iMikkola & ValtonenI [1991 lOuinlanI Il996al: iMerritL .200l1) . 
Unfortunately, neither p nor a are well defined for our A^- 
body models: p is formally divergent as r — > (eq. [TOl l. and o 
drops to zero a t the origin in the absence of the central binary 
(lDehnenll993l) . We can crudely evaluate equation ( fTSI l by set- 
ting p « 0.595(0.338) and a w 0.216(0.244), the mean and 
mass-weighted, rms values within a sphere of radius 0.1 (0.2) 
about the center of the (binary-free) Dehnen model. The re- 
sults, with H = 16, are s w 44(22). These are likely to be 
overestimates: the central density of the galaxy drops as the 
binary ejects stars and the central velocity dispersion is in- 
creased by the presence of the binary. If we decrease p by a 
factor of two to account for ejections and set o equal to the 
rms velocity dispersion in the y = 0.5 Dehnen model contain- 
ing a central, M = 0.01 point mass, the predicted hardening 
rates drop to ^ 13(8). These numbers are reasonably con- 
sistent with the low-A^ hardening rates shown in in Figure [8] 
s « 6.5, suggesting that the binary is approximately in the 
"full loss cone" regime at these low values of A^. 

4.2. Eccentricity Changes 

The A^-dependence of the eccentricity evolution (Figure |9]l 
is not quite so transparent. Although the two black holes 
were initially placed on circular trajectories, perturbations 
from passing stars sometimes resulted in very non-zero ec- 
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Fig. 7. — Binary hardening rate as a function of time, computed as an 
average over the independent integrations at each N. Line styles have the 
same meaning as in Figure |6] Tick marks indicate where the hardening rate 
was evaluated for Fig.[8] i.e., at {a-') =750. 

centricities developing around or even before the time the 
binary became hard. This was especially true in the small- 
A^ integrations (Fig. |5]); for A^ — ^K, the mass ratio between 
black hole and star was only 40 and a single star-binary in- 
teraction at early times could induce a substantial change in 
the binary's orbit. Once established at early times, these 
eccentricities tended to persist. Spurious changes in e in 
small-Af-body simulations h a ve been noted by other author s 
dOuinlan & HernquistI \l99% iMilosavlievic & Merritt! I200TI) . 
However the general trend in Figures |5] and |9] is clearly to- 
ward smaller eccentricities for larger A^. 

Statements about eccentricity evolution of massive binaries 
are often based on the results of three-body scattering ex- 
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Fig. 8. — W-dependence of the binary hardening rate, computed by fitting 
fl^'(f) in an interval A/ = 50 centered on a^' = 700. Asterices: W-body 
results (Fig. |6jb)), computed as averages over the set of ensembles at each 
N. Filled circles: Fokker-Planck results (Fig. I14t . omitting the "secondary 
slingshot." Open circles: Fokker-Planck results, including the "secondary 
slingshot." 

periments ("Mikkola & Valtonenll 991 lOuinla n'l 996a: ' Merritll 
12001 ). In these experiments, changes in e are typically ex- 
pressed in terms of changes in a as 



dt I \dt \a 



(17) 



where K — K{e,a) is a dimensionless rate coefficient and 
indicates averages over impact parameter and velocity at 
infinity. Mikkola & Valtonen (1992) and Quinlan (1996) 
give approximate analytic fits to Ki{e,a,Voo), the impact- 
parameter-averaged rate coefficient describing changes in e 
due to interaction of the binary with stars of a single velocity 
Voo. These expressions for Ki can be converted into expres- 
sions for K by averaging over an assumed velocity distribu- 
tion at infinity, and Quinlan (1996, Fig. 9) shows the results of 
such a calculation. (Sesana et al. 2006 present similar plots.) 
Evolution is always found to be in the direction of increasing 
eccentricity, i.e. K >0, excepting possibly in the case of soft, 
nearly-circular binaries (Quinlan 1996, Fig. 9d-f). Evolution 
rates tend to increase with increasing hardness of the binary, 
reaching maximum values of 0.2 for equal-mass binaries 
with e w 0.75 and falling to zero at e = and e = 1. This 
is at least qualitatively consistent with Figure |9] which shows 
de/dln{l /a) generally increasing at late times, i.e. for larger 
binding energies. 

Comparing these predictions quantitatively with the A^- 
body experiments is desirable, but problematic for a variety 
of reasons, the most important of which is probably the strong 
dependence of K on o/Vhi,,, where o is the stellar velocity dis- 
persion (assumed independent of position) and Vbi„ the binary 
orbital velocity. In galaxy models like ours, a is a steep func- 
tion of radius near the galaxy's center and it is not clear what 
value to choose. 

In the limit of large binding energy, Vbi„ ^ O, the velocity 
at infinity is irrelevant and K as determined by the scattering 
experiments becomes independent of a. Mikkola & Valtonen 
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Fig. 9. — Evolution of the mean value of e. Each fine is an average of the e 
values in the various A'-body integrations that started from different random 
realizatio ns o f the same initial conditions. Dashed lines show solutions to 
equation j20t . 

(1992) find for K in this limit the approximate expression 



(1 



2e 



(1-.2) 



2\ 



OT = 0.3e^-0.8 



(18a) 
(18b) 



while Quinlan (1996) gives, for an equal-mass binary in the 
large-binding-energy limit, 

K{e)Ke{l-e^)''° {ki+k2e), (19a) 
{hMM) = (0.731,0.265,0.230). (19b) 

Figure [To] shows that the two expressions are in good agree- 
ment. 

A rough value of o/Vhin in our simulations is ^ 2a^l^, 
where o has been set to ^ 0.2, its mean value within a 
sphere of radius 0. 1 (neglecting the effects of the binary). For 
fl^' in the range 500 < a^' < 2500 (Fig.|6ll, this expression 
gives 0.1 ^ ojVbin ^ 0.04. Figure 9 from Quinlan (1996) 
suggests that K{e) reaches its large-binding-energy limit for 
^I'^bin ^5 0.05, SO our simulations should be in or near this 
regime at late times. 

Accordingly, Figure|9]shows solutions to 



de = J^ie) dXna 



-1 



(20) 



using Quinlan's expression for K{e). The agreement with the 
A^-body results is quite reasonable, especially for the larger 
values of A^. Nevertheless, we stress again that the final ec- 
centricity values in our A^-body simulations are influenced 
strongly by noise-induced changes in e at early times, and 
these changes would be much smaller in the large-A^ regime 
of real galaxies. 

4.3. Mass Deficits 

As the binary hardens, it ejects stars from the nucleus and 
lowers its density. These density changes are sometimes es- 
timated from scattering experiments in a fixed background 
like those described above, e.g. the change in core mass 
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Fig. 10. — Two approximations, deri ved from three-body scattering ex- 
periments, for the coefficient K (eq. describing the rate of eccentricity 
evolution in the limit of large binding energy. Solid line: Quinlan (1996); 
dashed line: Mikkola & Valtonen (1992). The solid line was used to com- 
puted the evolutionary tracks (dashed lines) in Fig.|9] 

is equated with the mass "ejected" by the binary. How- 
ever the fact that the binary continues to harden at late times 
(Fig.|6j5) implies that depopulated orbits are continually being 
re-supplied. Changes in nuclear density are therefore a com- 
petition between ejection of stars (some of which may remain 
bound to the core) and re-population of orbits by gravitational 
scattering. A number of other mechanisms can also influence 
the evolution of the central density; for instance, loss of mat- 
ter from the core lowers its binding energy and causes it to 
expand. The net effect of these various processes is difficult 
to estimate without full A^-body simulations. 

We follow the standard practice of describing changes in 
core mass in terrns of t he mass deficit Mdef, defined by 
iMilosavlievic et an(l2002h as the difference in integrated mass 
between the density profile and the initial density profile, 
within the region influenced by the binary. Mass deficits 
have been estimated in a numbe r gala xies (Milosavlievic et al. 
12002^, 'Ravi ndranafli et al.ll2002l: [Graham 2004; Merritt 2006) 
using assumed forms for the pre-existing density profile. Fig- 
ure[TT]shows M^uf versus time, and versus binary semi-major 
axis, for the averaged A^-body integrations. 

As in previous work (IMilosavlievic & Merrittl20^ltlMerrit3 
12006'), we find that the mass deficit increases suddenly when 
a K, ah, to a value Mdef ~ M\2- Since the initial conditions 
adopted here are rather artificial - neither of the black hole 
particles was placed at the center, for instance - the value 
which we find for M^ef at this time may not accurately re- 
flect the value following a real galaxy merger. We therefore 
present in Figure \i\\Mri^f — M,uf i,. the change in the mass 
deficit since the time at which a = a/,; as above, we take this 
time to be f = 10 (Fig.|6ll. 

When plotted vs. a/Ja (Fig.fTTb). the A^-dependence of the 
evolution almost disappears, allowing the differential mass 
deficit to be expressed almost uniquely in terms of the change 
in semi-major axis. As shown below, a binary would not be 
expected to evolve past a^' w lOOa^^ before gravitational 
wave losses begin to dominate the evolution, implying a maxi- 
mum mass deficit of ^ SM^; however an extrapolation of this 
prediction to the much larger regime of real galaxies would 
be dangerous. Figure [12] shows averaged density profiles at 



various times for the integrations with A^ ~ 65K. 

5. THE FOKKER-PLANCK MODEL 

As shown in Figure [8] the A^-dependence of binary 
hardening rate in the A^-body simulations is s ^ N^^ . 
This is substantially flatter than the ~ A^^' dependence 
expected in a diffusively-repo pulated ("empty") loss cone 
(IMilosavlievic & Mertt3 120031) . which makes it difficult to 
extrapolate the A^-body results to the regime of real galaxies. 
In this section we develop a Fokker-Planck model that can re- 
produce the A^-body results and which can also be applied to 
systems with arbitrarily large A^. Unlike previous treatments 
of this problem based on encounter theory, we allow the radial 
distribution of matter to evolve in our Fokker-Planck models, 
due both to loss of stars that interact with the binary, and to 
diffusion in energy of non-interacting stars. These improve- 
ments will be shown to be crucial for accurately reproducing 
the A^-body results. They also allow us, for the first time, to 
make quantitative predictions about the evolution of the mass 
deficit in galaxies where binary evolution is driven by coUi- 
sional loss-cone repopulation. 

5.1. Loss-cone Dynamics 

Consider a spherical galaxy containing a massive central 
binary that acts like a sink, ejecting stars that come sufficiently 
close to it. Let E = — v^/2 + i|/(r) be the binding energy per 
unit mass of a star in the combined potential ^(r) = — ^(r) 
of the galaxy and the binary; the latter is approximated as 
—GMn/r. The binary defines a loss cone of orbits that satisfy 
J ^ Jic{E), where 



(21) 



here J is the angular momentum per unit mass of a star and 
ric is the radius of the ejection sphere around the binary. 

Suppose that the binary has interacted with and ejected all 
stars that were initially on orbits satisfying J < Ji^. (In Fig.[6l 
this appears to have occurred by a time of ~ 15.) The binary's 
subsequent hardening is limited by the rate at which stars are 
scattered onto previously depleted loss-cone orbits. A funda- 
mental quantity is the ratio qic{E) between the orbital period 
P{E) and the (orbit-averaged) time scale for diffusional refill- 
ing of the consumption zone (Paper I): 



qic{E) : 



1 



Ric{E) 



^lim«^ 

Vr R^Q 2R 



(22) 



Here R = J^/J.iEf is a dimensionless angular momentum 
variable, Q <R < 1, with Jc{E) the angular momentum of a 
circular orbit of energy £, and ((A/?) ) is the diffusion coef- 
ficient associated with R. The limit R ^ m equation ( |22] | 
reflects the approximation that only very eccentric orbits are 
scattered into the binary; the orbital period is likewise defined 
in terms of a 7 = orbit. This approximation breaks down 
for the most bound orbits but as we show, almost all of the 
loss cone repopulation comes from stars weakly bound to the 
binary. 

In the case of orbits with periods much shorter than the re- 
filling time iqic ^1), the system is "diffusive" and the loss 
cone is largely empty. For orbits with periods much longer 
than the refilling time {qic ^1), the system is in the "pin- 
hole" or "full loss cone" regime. In a galaxy containing a bi- 
nary with fixed r/^ , qic increases with decreasing E, i.e. with 
increasing distance from the binary. The energy at which 
qif. = 1 is defined as the critical energy, Ecrit, that separates 
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Fig. 1 1 . — Evolution of the mass deficits in tiie W-body integrations, vs. time (a) and semi-major axis (b). ai, is the binary separation at / = 10, when the hard 
binary forms, and Mjef.h is th^ mass deficit at this time. Line styles have the same meaning as in Figure|6] 




Fig. 12. — Evolution of the mean density profile in the 64K integrations. 
Black: / = 0; red: t = 50; blue: t = 150; green: t = 250; orange: t = 350. 

empty- from full loss cone regimes. The A^-dependence of the 
problem appears via the angular momentum diffusion coeffi- 
cient ((A/?^)), which scales (approximately) linearly with the 
mean stellar mass, i.e. inversely with for a fixed mass of the 
galaxy (Paper 1). Other factors that influence qic are the degree 
of central concentration of the galaxy (high central density 
implies larger qu) and the size r/^ of the interaction sphere 
(i.e. the binary semi-major axis). Milosavljevic & Merritt 
(2003) show that massive binaries in real galaxies (A^ ^ 10^) 
are essentially always in the empty loss cone regime, even in 
the extreme case of a p ~ r^^ stellar density cusp, due to the 
long relaxation times and to the large physical size of a binary. 

As a first step toward understanding the evolution of the bi- 
nary in our A^-body simulations, we plot in Figure [T3h qic{E) 
for our initial galaxy model, assuming two values for r/^. at 



each A^: rj^ = 100, corresponding to the time f w 15 when 
the hardening rate has just begun to exhibit a dependence on 
N (Fig.|6]l; and r/^ = a{t = 350), the final value of a (different 
for each A^). This figure suggests that none of the integrations 
was fully in the empty loss cone regime characteristic of real 
galaxies; even for A^ = 262k, qic > 1 except at energies close 
to \(/(r/,) (as defined above, rh is the gravitational influence ra- 
dius of the central mass, i.e. the radius containing a mass in 
stars equal to twice M12). As the binary hardens, qi^ increases 
in all of the simulations, and at the final time step, qi^ > 1 at 
E < x|/(r/,) for all A^, i.e. the binary has evolved essentially 
completely into the full loss cone regime. 

A more useful characterization of the binary's loss cone 
is shown in Figure \T3h . For this figure, the flux of stars 
into ric was computed, and broken into two parts: the flux 
Ffuii originating from stars at energies such that qi^. > 1 ; and 
Fempty, from Stars with energies such that qi^ < 1 . The energy- 
dependent flux !f {E) can be derived from the orbit-averaged 
equation describing diffusion in J (Eq. 19, Paper I); 



Ri, d fd9^\ 



where 9{,{E,R,t) ^ 4iz^P{E)J^{E)f{E ,R,t) is the number 
density of stars in the {E,R) plane. ^ The flux into the bi- 

' We assume in writing equation (23) that the orbit-averaged Fokker- 
Planck equation can be applied near the loss-cone boundary. This is valid 
for the diffusively-repopulated loss cone of a binary in a real galaxy, but may 
not be valid at low energies in the A'-body simulations since the loss cone is 
nearly full and the separation of time scales on which the orbit-averaging is 
based breaks down. Nevertheless equation j23) is traditionally applied even 
in this regime ( Cohn & Kulsrudtl978i ; ]MaEorrian & Tremainal999i) . Our ex- 
pression for the flux does tend to the correct limit in the full loss cone regime, 

2> 1 . See Shapiro & Marchant (1982) for a treatment of the loss cone that 
is not based on the orbit-averaged approximation 
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Fig. 13. — (a) The function cjidE) that describes the ratio of the orbital peiiod at E to the timescale for diffusional refilling of the loss cone; ^ S> 1 indicates 
that the loss cone is "full," and real galaxies have qic < 1. Line styles have the same meaning as in Figure|6] Thick curves show qic for o^' = 100, when binary 
has just entered the A'-dependent phase of its evolution (Fig. 6). Thin curves show qi^ for the binary at the final time step, t = 350; the binary separation at this 
time is different for each N. The radius of the loss sphere has been set to a. Vertical dotted line is £ = \(/(r;,). (b) The fraction of the flux of stars into the binary's 



loss cone that is contributed at energies where qic> 1, i.e., where the loss cone is essentially full. Lines show predictions for N = 
ignore binary-induced changes in the mass distribution of the galaxy. 



(0.5,1,2,4) X 10^ These plots 



nary is 



J {E)dE = 



d 
dt 

Ric 



9i{E,R,t)dR 



dE 



-qic 



R 



dH. 
dR 



dE 



4K^jUE)qi,{E) 



R 



dR 



dE. 



(24a) 
(24b) 
(24c) 



Ro 



In these expressions, / has been allowed to fall to zero at an 
angular momentum 7?o(£) that is different from /^/^.(ii). Cohn 
& Kulsrud (1979) derived an approximate expression for Rq: 



Rq{E)=R,,{E)x 



(exp{-qi^), qic{E) > 1 

\exp(-0.186^,,- 0.824^5), qi,{E) < 1. 



For qic ^ 1, /?o ~ Ric but as qic increases, the loss cone is 
largely full and Rq w 0. Finally, we adopt the steady-state 
solution to equation (|23] ) for /, i.e. 



-m 



in(i//?o)-r 

(assuming Rq <^l) implying a diffusive flux 

f{E) 



:r{E)dE^4n%.iE)q,,{E) 



ln(l//?o)-l 



dE. 



(25) 



(26) 



Here, / — Jq f{E ,R)dR is the isotropic / that has the same 
total number of stars at each E as the true f{E,R). 

As noted above, the loss cone of a binary black hole in a 
real galaxy is essentially empty, i.e. almost all of the stars 
scattered into the binary would come from energies E > Edt- 
Figure [T3b shows the results of applying equation ( |26] | to our 
initial A^-body model, with r;^ = a and with a allowed to vary 



over the range 100 < a < a 



-1 

r=350- 



In this figure, /y„// is 



the flux integrated from to Ec-u and F is the total flux. At 
the start of the A^-body integrations. Figure [T3b suggests that 
the binary in the larger-A^ models (A^ ^ 6Ak) is essentially in 
the empty loss cone regime, /y,,// ^ F . However by the final 
time step, the binary has shrunk and entered into the "pin- 
hole" regime, Ffuii > Fempty, for all A^. In the integrations 
with A^ ^ 16^, the binary is in the full loss cone regime from 
the start. 

Figure [Tjb also includes curves for the cases A' = 
(0.5,1, 2, 4)x 10^. Values of A^ up to 4 X lO*" are now compu- 
tationally feasible via direct-su mmation code s combined with 
special-purpose hardware (Har fst et al.ll2006h . and Figure[T3b 
suggests that this A^ value is large enough to place the binary 
effectively in the empty loss cone regime for most of its evo- 
lution. (The minimum required A^ would be larger than this if 
the binary were given the smaller mass, ^ lO^^Mg^i, typical 
of black holes in real galaxies, or if the galaxy model were 
more centrally concentrated.) 

Figure[T3]illustrates the difficulty of scaling the binary evo- 
lution observed in our A^-body simulations to real galaxies. In 
the empty loss cone (i.e. diffusive, large-A^) limit, the sup- 
ply of stars to the binary scales as ( (A/?^) ) <:x: °^ A^-i for 
a fixed total galaxy mass (ignoring the weak dependence of 
the Coulomb logarithm on A^). In the full loss cone (pin- 
hole, small-A^) limit, the loss cone flux is independent of A^. 
In between these limiting cases, one expects (Paper 1) that 
the flux, and hence the hardening rate of the binary, scales as 

A^^P, < p < 1. This is consistent with the s ^ A^-o.36 
pendence observed here (equation [ISTl. Figure [T3b suggests 
that of order A^ w 10^ stars would be required before the bi- 
nary is comfortably in the empty loss cone regime, allowing 
its evolution to be reliably scaled to larger values of A^. 

Even if we were in this regime, the expressions given above 
for the flux of stars into the binary's loss cone might not 
accurately predict the binary's evolution, since they ignore 
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changes in the galaxy's structure. Figure[T2]suggest that these 
changes are significant: the density near the galaxy's center 
changes by a factor ~ 2 as the binary hardens. We now con- 
sider a model that includes both changes in the binary due to 
interaction with stars, as well as binary-induced changes in 
the stellar distribution, and that can be reliably scaled to the 
large-A^ regime of real galaxies. 

5.2. Evolutionary Model 

The 7-directed flux of stars into the binary, described by 
equation ( |26] |. implies a decrease in the number of stars at 
Jic Jc{E). In Paper 1, this decrease was followed by in- 

tegrating equation ( |23] ) forward in time at fixed E. The justi- 
fication for treating the problem in this restricted way was the 
difference in time scales between E- and /-diffusion; the for- 
mer occurs in a time ^ TJ- while the latter requires ~ {ai,/r)Tj. 
The evolution of f{J;E) and f (E) over the shorter of these 
time scales was followed starting from a completely emptied 
loss cone and the change in the density of the core was com- 
puted from the changes in N{J;E) at every E. 

In the present paper, we focus on changes that take place 
over the longer of these two time scales, ^ Ti-. This allows 
us to largely ignore the initial conditions, and to assume that 
an expression like (l25T l is an adequate description of the J- 
dependence of / at every E. However it also implies that we 
can not ignore changes in E, which occur on timescales of 

On these longer time scales, the evolution of the density 
near the binary is a competition between loss of stars that dif- 
fuse onto low-7 orbits and are ejected by the binary, as de- 
scribed by J (£■), and replenishment due to stars that diffuse 
in energy from regions of lower E, i.e. larger radius. Be- 
yond a certain radius, the relaxation time is so long that the 
/i-directed flux can not compensate for the integrated loss- 
cone flux, / f {E)dE, and the mean density within this radius 
must drop - implying the creation of a mass deficit. 

We can approximate the evolution of the galaxy/binary sys- 
tem in this late-time regime via a modification of the orbit- 
averaged Fokker- Planck equation for f{E): 

dN dpE , , 



a? 'dE 

where J is the /-directed flux defined in equation 
Fe is the energy-directed flux, given by 



Fe 



Dee 



3/ 



:647i%2m*lnA 



dE'q{E')f{E') 



D£ = -647i'^G^m*lnA 



(27) 
, and 

(28a) 



(28b) 

' dE'p{E')f{E') (28c) 



q{E) dE'f{E') 



In these expressions, f{E) is understood to be the mass 
density of stars in phase space associated with the function 
/(£) defined above, and the quantities Fe and J are mass 
fluxes. N{E)dE ^4K^p{E)f{E)dE is energy-space distribu- 
tion, with p{E) and q{E) the phase-space weighting factors. 



p(£)=4 



AE) 







v{r)r^dr, 



3 JO 



AE) 



v^{r)r^dr, 



(29a) 
(29b) 



1 /2 

and V = [2<t>(r) — 2E] ' . Near the binary, where the poten- 
tial is close to Keplerian, p{E) « 2-^IHG^M\2\E\-^I^ and 
q{E) = (2'/27i/6)G^M^2|£|-3/2. InA w ln(Mi2/mJ is the 
Coulomb logarithm. Henceforth / and are explicitly de- 
fined as mass (not number) densities, and J is the mass flux 
into the binary's loss cone. 

An equation like ( |27] |. in which the 7-dependence of / 
is contained implicitly in J {E,t), was first written by Bah- 
call & Wolf (1977). It has since been adopted by a num- 
ber of other authors to describe the evolution of the distribu- 
tion of stars, compact obj ects or dark matter around a single 
supermassive black hole dMurphv et all 1 1991; Merritt 2004 
iHopman & Alexander ''2Q0o). It is being used for the first time 
in the present paper to describe the evolution of the stellar 
distribution about a binary black hole. Since f scales only 
as ^ log/"/^.' (equation |26]|. the ratio of the two terms on the 
right hand side of equation ( |27] i is not greatly affected by the 
much greater (in linear extent) size of the loss cone of a binary 
compared with a single black hole. 

The relation between the flux into the binary's loss cone and 
the rate of change of its semi-major axis a is 

with /J = M1M2/M, the binary reduced mass and AE{E) the 
mean specific energy change of stars, originally at energy E, 
that interact with the binary. In Paper 11, we set 



AE = AE, 



Hills 



-(c) 



(31a) 



s{t) = — - 
dt \a 



m 

aM 



I !F{E,t)dE. (31b) 



The coefficient (C) is independent of energy for stars that in- 
teract wath_aj'h^d^ (Hills 1983; Mikkola & Valtone^ 
1 1 9921; l0uinlan|[T996al) . Hardness is defined as Vft,„/a, where 

Vbin = \/ GM 12 /a is the relative velocity of the components of 
the binary and a is the stellar velocity dispersion in the unper- 
turbed galaxy. A n equal-m ass binary is in the "hard" regime 
when Vbin/o ^ 3 (IOuinlantl996a) . In the current models. 



Vbi, 

On 



■0.36 



1/2 



(32) 



with Op ~ 0.278 the peak velocity dispersion in the y = 0.5 
Dehnen model. The A^-dependent phase of binary evolution 
begins at a^' « 100 in the A^-body models (Fig. |6]l, hence 
Vbin/Op ^ 3.6 and equation ( l3Tl i is expected to be accurate. 
In Paper 11, setting (C) w 1 .25 was found to reproduce the A^- 
body hardening rates. Yu (2002) argued for a similar value of 

Expressions like ( 131 al l were derived from scattering experi- 
ments that allowed for the possibility of multiple interactions 
between star and binary. However the confining effect of 
the galaxy's gravitational poten tial was ignored. As noted in 
iMilosavlievic & Merritt! (l2003b . stars ejected once by the bi- 
nary can interact with it again as they return to the nucleus on 
nearly-radial orbits. If the energy change during the first in- 
teraction is not large enough to eject the star completely from 
the galaxy, it will experience one or more "secondary sling- 
shots", and the total energy extracted from the binary by the 
star will be the sum of the discrete energy changes during the 
interactions. 
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Fig. 14. — (a) Fokker-Planck evolution of binary semi-major axis in a set of integrations designed to mimic the A'-body simulations with N = 65k. Dot-dashed 
line: fixed potential and density; dashed line: fixed potential, evolving density; thin solid line: evolving density and potential, no re-ejections; thick solid line: 
evolving density and potential with re-ejections, (b) Fokker-Planck integrations with parameters chosen to mimic the W-body simulations with various N; color 
coding is the same as in the A^-body figures above. Thin lines: no re-ejections; solid lines: with re-ejections. 
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Fig. 15. — Evolution of the mass deficit in the suite of Fokker-Planck integrations presented in Fig. 1 141 Line styles have the same meaning as in Figure|6] 



A minimum condition for re-ejection is that a star remain 
bound to the galaxy after its first interaction with the bi- 
nary, E + AE ^ 0. Most s tars that interact with the bi - 
nary have apocenters ^ r/, ( Milosavljevic & Merrittll2"003h : 
since the gravitational potential at this radius is dominated 
by the galaxy, we can write this condition for re-ejection as 
\AE\ <I>j,(0) with 4>*(0) the central value of the galaxy's 
(stellar) gravitational potential. The y = 0.5 Dehnen models 
used here have 't>*(0) « 0.67 in the adopted units, implying 
a"' ^ 200 for re-ejection. 

Even if a star satisfies this condition, re-ejection will only 



be effective if the star remains in the binary's loss cone for 
longer than an orbital period, i.e. if q{E) <S 1. Re-ejection will 
also fail for a star with apocenter greater than some r,„„i S> o,, 
since the overall potential in a real galaxy is never precisely 
spherical and the star will be perturbed frorn its nearly radial 
orbit on the way in or out Vicari et al] (l2006h . 

We considered a modified form of equation (ISTT l that 
accounts for re-ejections. Let AZsmax = 't'('"mo.x) — '^{fh)- 
Re-ejection was assumed to occur if the following condi- 
tions were both satisfied: (i) |AZiHiiis| < [A/ima^l; (ii) q{E + 
AZiHiiis) < Qmax ~ 1- Condition (i) guarantees that the star 
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remains bound to the galaxy after the first ejection, with 
apocenter < rnuu- This condition is roughly equivalent to 
Gfi/a < 4>ga/(0) « 1 and is satisfied for a"' <, 500 « 5aJ^^ 
in our models, i.e. during the early phases of binary evolu- 
tion. Condition (ii) guarantees that the star will remain within 
the binary loss cone for of order one orbital period or longer 
after the first ejection; this condition is satisfied at large (i.e. 
bound) values of E (Fig. 13). 

For fl^' greater than ^ 160, which occurs shortly after for- 
mation of a hard binary (Fig. 5), even a single re-ejection 
would give a star enough energy to escape the galaxy. Hence, 
at energies such that conditions (i) and (ii) were both satis- 
fied, we set AZi = 2A£'Hiiisi while if either condition was not 
satisfied, re-ejection was assumed not to occur and we set 
AZi — A/iHiiis- This scheme has two parameters, q,„ax and r„^ax', 
the results are weakly dependent on r^ax for r^ux ^ Oi and we 
fixed r,nax — lOOr/,. The consequences of varying q,„cix are dis- 
cussed below. 

Finally, we need to account for changes in the gravita- 
tional potential as the stellar distribution evolves. (We ig- 
nore possible changes in the mass of the binary.) Here we 
follow Henon's (1961) scheme of assuming that / remains a 
fixed function of the radial adiabatic invariant as the poten- 
tial is adjusted. Our numerical schemes for advancing f was 
based closely on the algorithms described bv lCohnI (119801) and 
IOuinlanl(ll996bl) . 



galaxy model was a Dehnen sphere, equation ( fTOb . with y = 
0.5. This is the shallowest central slope that is consistent 
with an isotropic phase-space distribution around a central 
point mass; it is also a fair representation of the core pro- 
files that are produced during the "rapid" phase of cusp de- 
struction tha t accompanies the ini tial formation of the mas- 
sive binary dMerritt & Szelll l2006l) . Fokker-Planck integra- 
tions were carried out for different values of N = Moai/m^, ~ 
(lO'', 10^, 10^^). The time axis in these plots is the relax- 
ation time measured at the binary's influence radius in the ini- 
tial model; all integrations were continued until t — 47; (r/,). 
Equation ( fT2] ) was used to set the initial value of a; quantities 
like the mass deficit in Figures[T6]and[T7]should be interpreted 
as the accumulated change in these quantities after the binary 
first becomes "hard." Unless otherwise stated, re-ejections 
were ignored. 

In all of the integrations, the binary begins in the diffu- 
sive, or empty loss cone, regime {qu- » 1) due to its large 
initial separation, and evolves toward the pinhole, or full loss 
cone, regime {qi^ > 1) as it hardens. The transition to the 
pinhole regime occurs later for larger A^; for = lO'^ (the 
heavy curves in Figs. [TSlandfTTb the binary remains essen- 
tially in the diffusive regime until the end of the integration at 
4Tr{ri,). However we argue below that evolution of binaries in 
real galaxies would typically be expected to terminate before 
the pinhole regime is reached. 



5.3. Comparison with the N-Body Integrations 

Figure [TJa) compares the evolution of a^^ in a set of 
Fokker-Planck integrations with initial conditions chosen to 
mimic those in the A^-body integrations (y = 0.5, M ~ 0.01, 
fl"' (f = 0) = 0.01). Fixing p(r) and <I>(r) (dot-dashed curve) 
is equivalent to the assumptions made by Yu (2002), who ig- 
nored changes in the stellar distribution as the binary evolved. 
Allowing the density and potential to evolve (solid lines) re- 
sults in a considerably lower hardening rate for the binary. 
Including the secondary-slingshot (heavy solid line) increases 
the hardening rate but only slightly; as explained above, once 
the binary becomes hard, most stars that interact with it are 
ejected completely from the galaxy and do not return to the 
binary's sphere of influence. 

Figure fT4T b). which can be compared with Figure |6jb), 
shows the evolution of binary semi-major axis in a set of 
Fokker-Planck integrations with the same values of as in 
the A^-body integrations. The correspondence is quite good; 
the Fokker-Planck integrations show a slightly steeper depen- 
dence of the binary hardening rate on A^ (Fig. |8]l. The evo- 
lution of the mass deficit as derived from the Fokker-Planck 
integrations is shown in Figure [15] (cf. Fig. fTTT i. Here the 
correspondence is not quite as good, but still reasonable; the 
weak dependence of Mdef on A^ for large A^ is well reproduced. 

6. PREDICTIONS OF THE FOKKER-PLANCK MODEL FOR LARGE N 

Having established that the Fokker-Planck model can 
mimic the joint binary/galaxy evolution seen in the A^-body 
integrations, for various values of N <, 10^, we now extend 
this model to the much larger A^ regime of real galaxies. The 
goal is both to predict the long-term evolution of a massive 
binary in a real galaxy, and also to record the changes in the 
central structure of the galaxy. 

Results for a galaxy containing a binary with M12 = Mi + 
M2 = lO^^'Mgai and two mass ratios a = M2/M1 = (1,0.1) 
are shown in Figures [16] and [17] respectively. The initial 
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6.1. Binary Hardening Rates 
]and[T7]show that at large A^, the binary hardening 



lard ■ 



(33) 



tends to a fixed fraction of Tr{ri,) at any given a. This is 
the "empty loss cone" regime. The ratio 7hard/7i ('"/!) in- 
creases from 0.1 at large a, i.e. early times, to ^ 0.3 when 
a lO^^fl/,, with a weak dependence on binary mass ratio. 
We will argue below that binary black holes in real galaxies 
lie close to the large-A^ hardening curves throughout much of 
their evolution and so it is of interest to develop an analytic 
understanding of this regime. 

Since qic{E) (equation |22]) is the ratio of the orbital pe- 
riod to the diffusional loss cone refilling time at energy E, 
i.e. qic{E) w P{E) /[Ric{E)Tr{E)], we can rewrite the flux of 
stars into the binary, equation (|26] |. as 



T {E)dE w 4k^J^{E)P{E)T-'{E) 



-I, 



f{E) 



ln(l//?o)-l 



dE. (34) 



Assuming a fixed mass model for the galaxy, the flux into the 
binary, integrated over one relaxation time, scales therefore as 



1-1 



r{E)T,{E)^[\n{\/RQ)-\^ 



(35a) 
(35b) 
(35c) 



The binary hardening rate is fixed by J and a (equation l31bl ). 
so these expressions imply that the binary's evolution over 
a specified number of relaxation times will be smaller for 
smaller A^, i.e. larger qu ', while in the large-A^ limit, the evo- 
lution rate at a given a will be determined solely by T^. These 
predictions are consistent with the upper panels of Figures [T6l 
and[T7] 
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Fig. 16. — Joint binary-galaxy evolution in Fokker-Planck models with Mi = M2 and M12 = lO^^Mg-^i. (a) Binary semi-major axis; (b) binary hardening time; 
(c) mass deficit as a function of time, and (d) mass deficit as a function of binary separation. Different lines correspond to different values of N = Mjai/m^t: 
A? = 10*, 10', 10" , 10'^ (thick line). Symbols mark the time t^^ at which the binary hardening rate equals the gravitational radiation evolution rate, assuming a 
binary mass of 10^ (squares), 10* (circles), 10^ M,., (triangles), and lO^M ;, (stars). Filled symbols denote models in which N is roughly equal to its value 
in real galaxies, for each value of M, . Dashed lines in panels (a) and (b) are the analytic model described in the text. 

The gradual decrease with time of the hardening rate is due 
to two factors: the decreasing size of the binary, and the de- 
cHning density of the core. Again ignoring changes in the 
core structure, the expressions given above can be used to es- 
timate how the hardening rate varies with a. The result, in the 
large-A^ limit, is 



for a = (1,0.1) respectively. The weak dependence of the 
fitting parameters on binary mass ratio reflects the lack of a 
mass ratio dependence in the evolution equations OTT l. 

Integrating equation ( |37] | gives a simple expression for the 
time dependence of the binary semi-major axis: 



^ard 



In 



(36) 



In 



(?) 




A Tr{rh) 



(39) 



i.e. the fractional change in a over one relaxation time is 
weakly dependent on a for large A^. 

We tried fitting a similar function to the large-A^ hardening 
curves in Figures [16] and [17] i.e. 



1 



= Aln(^ 



B. 



(37) 



The results are shown as the dashed lines in Figures [T6b 
and[T7b. We found good fits for 



A = (0.016,0.017), B= (0.08,0.09) 



(38) 



where t is defined, as in Figures [16] and [T7] as the time since 
the binary first became hard, i.e. the time since a — a/,. This 
function is plotted in Figures [T6h and[T7h. where it again pro- 
vides an excellent fit to the large-A^ evolution curves. 

We now show that real black hole binaries are expected to 
be in this empty loss cone regime throughout most or all of 
their evolution. Maximum traversal of the tracks in Figures[T6l 
and [TT] will occur if no physical process, aside from interac- 
tions with stars, affects the hardening rate until the gravita- 
tional radiation regime is reached. The time scale associated 
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Fig. 17.— LikeFig.[T6]butforM2 =0.1Mi andM|2 = IQ-^Mg^i. 

with gravitational radiation is (iPetersll 19641) 

a gr u-r w /JIVI j2 

where /j = M1M2/M12 is the reduced mass of the bi- 
nary and a circular orbit ha s been assumed. Following 
iMerritt & Milosavlievid (l2005b . Tgr can be expressed in terms 
of M, = M12 and ah using the M, — a relation, equation (O, 

as 

7g,«5.7xlO>^^M.-;f5al2 (41) 

with a = a/ai, and fl_2 = a/(0.01fl/,). 

We define teq as the time when T^ard = 7gr- In order to ex- 
tract teg in physical units from the Fokker-Planck integrations, 
we need to assign a value in years to 7; (r/,). This we do via 
the straight-line fit to the data in Figure [T] Combining equa- 
tions (IHi, (|37] |. and (HTt the condition Thaid = Tgi- becomes 

, 4 



= 7.1 X ioV(i 



2.19 



(42) 

with M,.6 = M,/1O^M0. We define aeq as the value of a that 
satisfies this equation. ForM. (10^, 10'', 10^, 1O^)M0, i.e. 



ogio (Oh/a) 



a « (44,70,112,180) km s-1, and using the values of A and 
B derived above, equation (l42T i implies 

fl/,/«eq« (315,93,27,8.0) (43) 
for a = M2/M1 — 1, and 

a,,/fleq~ (140,40, 12,3.5) (44) 

for a = 0.1. The corresponding times are 

(0.73,0.53,0.35,0.20) x 7;(r/,)(a = 1), (45) 
w (0.65,0.54,0.27,0.13) x 7;.(r/,)(a = 0.1). (46) 

The values just computed for Oeq and t^g correspond to 
the large-A^ (empty loss cone) limit of the Fokker-Planck 
equation, i.e. to the heavy curves in Figures [16] and [17] 
The filled symbols in those figures show where Thaid = Tgr 
on the four tracks that best correspond to the four val- 
ues just considered for M,. Since M. sa 1 x IQ^^Mgai 
(IMerritt & FerraresdlMTl) . we set N = (10^ 10^ 10^°, lO") 
for M. = (10^ 10^ 10^ 1O^)M0. The symbols confirm that 
binary black holes of mass M, > IO^'^Mq remain essentially 
in the empty loss cone regime throughout their evolution. For 
binaries of mass M, = IO^Mq, the evolution just prior to the 
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gravitational radiation regime begins to depart from that of 
a diffusive loss cone, resulting in somewhat lower hardening 
rates than predicted by equations (|37] | and ( [39] l. The discrep- 
ancy with the analytic expressions would be expected to in- 
crease still more for binaries of still lower mass (if such exist). 

6.2. Mass Deficits 

Next we consider the effect of the binary on the structure 
of the galaxy's core. Evolution of mass deficits is plotted in 
the lower panels of Figures [16] and [17] Particularly striking 
are the panels showing M^gf vs. binary hardness, fl/,/fl. As 
was true in the A^-body integrations, long-term evolution of 
the binary generates mass deficits that are very well predicted 
by the change in binding energy of the binary black hole, i.e. 
by fl/,/fl. This dependence is accurately described by 



Mi2 



(1.8,1.6)logio (ah/a) 



(47) 



where the numbers in parentheses refer to a ( 1 , 0. 1 ) respec- 
tively. The mass deficits generated between formation of a 
hard binary, and the start of the gravitational radiation regime, 
are given by setting a = aeq in this expression, i.e. 

Mdef« (4.5,3.5,2.6, 1.6)Mi2 (a=l) (48) 
w (3.4,2.6, 1.7, 0.9)Mi2 (a = 0.1). (49) 

for Mn = (10^ 10^ 10^ 10^)Mq. These values should be 
added to the mass deficits Mdef.h generated during the rapid 
phase of binary formation, i.e. Mdefh ~ 0.7a°-^Mi2 (iMerritj 
l2006l) . 

Mass deficits in these models are not related in a simple way 
to the mass in stars "ejected" by the binary. The flux of stars 
into the binary constitutes a loss term, — f {E,t), on the right 
hand side of equation ( |27] |. and in the absence of any other 
influences, the density of stars near the center of the galaxy 
would drop in response to this term. Removal of stars also 
reduces the gravitational force near the center, contributing to 
the expansion. However the second term on the right hand 
side of equation (i27] i. —dpE/dE, has the opposite effect. This 
term represents the change in N{E,t) due to diffusion of stars 
in energy; as the mass deficit increases, so do the gradients in 
/, which tend to increase the energy flux and counteract the 
drop in density. 

In principle, these two terms could balance, at least over 
some range in energies, allowing the binary to harden without 
generating a mass deficit. This would require 



FEiE) 



7 {E)dE, 



(50) 



i.e. the inward flux of stars due to energy diffusion at energy 
E must equal the integrated loss to the binary at all energies 
greater than E. However, at sufficiently great distances from 
the binary, the relaxation time is so long that the local Fe{E) 
must drop below the integrated loss term, implying that the 
density within this radius will drop. Growth of a mass deficit 
reflects the imbalance between these two terms. 

We illustrate this imbalance in Figure [18] which shows 
Fe{E) and / f {E)dE in the Fokker-Planck integration with 
a = 1 at a time ~ 7J (r/,). The lowest energy in the figure 
corresponds roughly to the outer edge of the binary-generated 
core. 

Yet another mechanism contributes to the growth of mass 
deficits in the Fokker-Planck models. Even in the absence of 
the loss term associated with the binary, the nuclear density 




0.59 0.6 0.61 0.62 0.63 0.64 0.65 



Fig. 18. — Fluxes in the Fokker-Planck integration withA'= 10 anda= 1, 
at a time ~ Tr{ri,). 

profile adopted here for the initial models, p ~ r^*^^, implies 
a "temperature inversion," i.e. a velocity dispersion that in- 
creases with radius. Relaxation drives such a nucleus toward 
a locally "isothermal" form before the onset of core collapse, 
causing the central density to drop (Ouinlan 1996b). The bi- 
nary contributes to this process by maintaining a flat density 
profile near the center, forcing the temperature inversion to 
persist. 

7. IMPLICATIONS FOR BINARY EVOLUTION IN GALAXIES 

Equation ( |39l ), based on the Fokker-Planck integrations, 
accurately describes the evolution of a hard binary in the 
empty loss cone regime (i.e. in galaxies with Mn ^ lO^'^M©) 
given the relaxation time at r/,, while equation ((4]i, based 
on observed properties of galactic nuclei (Fig. [TJ, provides 
the mean value of Tr{ri,) for galaxies with black hole mass 
M, ^ M\2- Together with equation ( i40b for the gravitational 
radiation time scale, these relations can be used to predict 
mean evolution rates and binary separations in real galaxies 
given (M.,a) = (M2/M1). 

Including the effect of energy lost to gravitational radiation, 
the binary's semi-major axis evolves as 



d n\ _ d n 

dt \ a J dt \ a 



hard 



- - 

dt \ a 



-T^\a) 



The time for the separation to drop from a^ to a is 



where 



10' Vx 



C=1.25M, 



Ay + B 



C + D{Ay + B)e'^y 



-1.54 



£)=1.75 X 10"'^a"^(l 



(51) 
(52) 

(53) 

(54a) 
(54b) 



andy,„ax ~ ln(fl/,/fl). The full time to coalescence, fcoab start- 
ing from fl/, is given by setting y,„ax = °o in this expression. 
Figure [19] shows t^oai as a function of M\2 for a = (1,0.1). 
Shown separately on this figure is the time to evolve from 
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(Mq) 

Fig. 19. — Time to coalescence starting from a = ai, as a function of 
binary mass. Solid curves are derived from equation i53\ with 
black/thick: a = 1 ; blue/thin: a = 0. 1 . Dotted curves show the evolution time 
from a = Oeq to a = 0, i.e. the time spent in the gravitational radiation regime 
only. Equation |55| gives accurate analytic approximations to fcoai{Wi2;0'). 

fl = fleq to fl = 0, i.e. the time spent in the gravitational ra- 
diation regime alone. The latter time is a factor ~ 10 shorter 
than the total evolution time fcoai, which motivates fitting the 
following functional form to koaiiMn'M): 



Y = Ci+C2X + C3X'-, 

^coal 



F = log 



10 



IQiOyr 



(55a) 
(55b) 

(55c) 



(This functional form is the integral of equation |39]) A least- 
squares fit to the curves in Figure[T9l gives 

a = 1 : Ci = -0.372, C2 = 1.384, C3 = -0.025 (56a) 
a^O.l : Ci = -0.478, d = 1.357, C3 = -0.041 (56b) 

The fit of the analytic expressions is better than 2% (a = 1) 
and 5% (a = 0.1); most of the deviations occur at the high- 
Mi 2 end where coalescence times are much longer than a 
Hubble time. 

Based on Figure [T9l binary black holes would be expected 
to reach gravitational wave coalescence in 10 Gyr in galaxies 
withMi2 ;$2 X lO'^M©. 

Figure |20] shows the probability predicted by equation ( |52| | 
of finding the binary in a unit interval of In a, 



P(lnfl) oc a 



da 
dt 



T{a), 



(57) 



for four values of M12 and for a = (1,0.1). Viewed at a ran- 
dom time before coalescence, a hard binary is most likely to 
be seen at a w 2fleq, although the distributions are nearly flat 
for 1 < aij/a <^ lah/agq. ForM^ ^ lO^M©, evolution for 10 
Gyr would only bring the binary separation slightly below o/,; 
in these galaxies the most l ikely separati on to find a binary 
would be the stalling radius (lMerrittll2006h . 
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Fig. 20. — Probability of finding a binary black hole in a unit interval 
of Ina. From left to right, curves are for Mn = (0.1, 1, 10, 100) x lO^M,;,. 
SoHd(dashed) curves are for M2/M1 = a = 1(0.1). Open circles indicate 
; filled circles correspond to an elapsed time since a = ai, of lO'" yr 



For the two smallest values of M, 
light. 



, the latter time occurs off the graph to the 



8. IMPLICATIONS FOR THE STRUCTURE OF GALAXY CORES 

In the most luminous spheroids, mass deficits generated by 
a binary black hole are likely to persist for the lifetime of 
the galaxy, since relaxation times are much too long for star- 
star scattering to alter the phase-space density (cf. Fig. [T]). 
In collisional nuclei on the other hand, relaxation times are 
short enough that the stellar distribution can be substantially 
affected by gravitational encounters after the binary black 
hole has coalesced into a single black hole. A Bahcall-Wolf 
(1976) cusp will form in a time ~ 7'r('"/i) after the binary 
black hole coalesces in to a single hole, inside a radius ^ 0.2r/, 
jMerritt & Szelll 12006 ). In addition, the structure of the nu- 
cleus beyond the cusp will continue to evolve, as two-body en- 
counters drive the stellar "temperature" profile toward isother- 
mality prior to the onset of core collapse (Ouinlan 1996b3. 
The nuclear density profile at some time after coalescence will 
depend on how far along the evolutionary tracks of Figs. [1§] 
and[T7]the binary evolved before coalescing, as well as on the 
elapsed time since coalescence. 

Figure |2T| illustrates these competing effects with a con- 
crete example. A Fokker-Planck integration with = 10^ and 
a = 1 was carried out until a time f = teq', teq was computed 
as above assuming a binary mass of lO^M©. The binary was 
assumed to become a single black hole at this time; the in- 
tegration was then continued for a time Ti-{rf,), but with the 
binary loss term f {'E) set to zero. As the figure shows, a 
p ~ r^^/^ cusp is generated at r < 0.2r/,. The net result is a 
flat core containing at its center a compact star cluster around 
the black hole. The stellar mass within the cusp is ~ O.IM,. 

If binary coalescence were assumed to take place sooner 
than ^ teq (due e.g. to gas-dynamical torques), the mass 
deficit would be smaller than the value ^ 3.5Mi2 generated 
in this integrati on, resulting in a nu clear density profile more 
like those of IMerritt & Szelll ([2006). As shown in that pa- 
per, a regenerated cusp can closely approximate the (coreless) 
density profile at the center of the Milky Way if the elapsed 
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Fig. 2 1 . — Stellar density and mass profiles in a Fokker-Planck integration 
with N = lO'. The binary black hole was assumed to coalesce at ( = 
(based on an assumed binary mass of lO^M,;,) and the integration was then 
continued, without the binary sink term, for one relaxation time at r;,. A 
Bahcall-Wolf cusp is generated at r < 0.2r;, ; the stellar mass within the cusp 
is ~ 0. IM. . Dashed line is the initial galaxy model. 

time since binary coalescence is ^ 8 Gyr In the integration 
of Figure |2T1 on the other hand, the larger mass deficit is not 
completely "erased" by formation of the cusp. 

A nuclear cusp like that in Figure |2T] would be unre- 
solved in all but the nearest galaxies. In fact, recent obser- 
vations suggest the presence of compact stellar nuclei ("nu- 
clear star cluste rs") at the centers o f most spheroids fainter 
than ~ lO'^Lr^ (iRossa et alJ l2006t IWehner & Hrnisl 120061: 
iFerrarese et al.ll2006h . The mean mass associated with the nu- 
clei is a fraction ~ 0.2% that of the host gala xy with a ±la 
range of 0.06% -0.52% (IFerrarese et al.ll2()06l) . If we assume 
that low-luminosity spheroids contain massive black holes 
and that the ratio of black hole mass to spheroid mass is sim- 
ilar to the mean value ^ 0.12% charact eristic of more lumi- 
nous galaxies ( Merritt & Ferraresel2001l) . the observed nuclei 
would have masses that are fractions 0.5 — 4 that of the black 
holes. This is somewhat larger than the value Mcusp/M, « 0. 1 
in the example of Figure |2Tl on the other hand it is possible 
that black holes in faint spheroids carry a larger fraction of the 
spheroid mass. The compact nuclei might also form in very 
different ways, e.g. from gas that accumulates at the center. 

9. ALTERNATE MODELS FOR BINARY EVOLUTION 

rv\il (120021) computed evolutionary tracks for binary black 
holes at the centers of a sample of early-type galaxies for 
which detailed luminosity profiles were available. Evolution 
beyond a k, ai, was modelled using the second term on the 
right hand side of equation ( |27] i. i.e. the /-directed flux of 
stars into the binary. The stellar distribution function was as- 
sumed fixed in time; binary-induced changes in the structure 
of the nucleus were ignored, as were changes in stellar energy, 
although the time scales associated with both sorts of change 
are comparable to the time scale for loss cone repopulation. 
As shown here (FigureO, allowing for changes in the struc- 
ture of the nucleus in a Fokker-Planck model reduces the bi- 
nary hardening rate by a factor ~ 2. Yu (2002) concluded that 
binary black holes in spherical galaxies with o < 90 km s^^ 
could coalesce in a Hubble time. This velocity dispersion cor- 
responds to a binary mass of- 3.5 X IO'^Mq (eq.O. Yu's con- 



clusion is consistent with, but slightly more optimistic than, 
the one reached in the current study (see Fig.[T9b: the differ- 
ences are probably due to Yu's neglect of the back-reaction of 
the binary on the nucleus. As in the current study, Yu found a 
weak dependence of coalescence time on binary mass ratio. 

Some recent studies have inferred rapid evolution of super- 
massive binary black holes at the centers of spherical galax- 
ies, even in the absence of collisional loss-cone repopula- 
tion. ISesana et alJ (l2007h used detailed three-body scattering 
experiment s to evaluate th e effectiveness of the "secondary 
slingshot" (Milosav lievic & Merritt 2003) at extracting en- 
ergy from massive binaries after they had reached the stalling 
radius a « a/, in spherical galaxies. They found that binaries 
could shrink beyond a/, by factors of — 4(2) for mass ratios 
of 1(0.1); for mass ratios below ^ 0.01 the secondary sling- 
shot was found to be ineffective. Almost all of this evolution 
took place within a few galaxy crossing times after the hard 
binary had formed; after this time, all of the stars that were 
originally within the binary's loss cone had been completely 
ejected from the galaxy. 

In spite of this very modest evolution, I Sesana et al.l ( l2007h 
concluded that "even in the absence of other mechanisms 
driving orbital decay, pairs involving genuinely supermassive 
holes [i.e. with combined mass ^ lO^MpJ should not stall". 
This optimistic conclusion appears to have been based on an 
evaluation of the mass "ejected" by the binary (their Fig. 5), 
rather than on the more fundamental criterion of binary sep- 
aration. The time to coalescence once a binary reaches the 
gravitational radiation regime is 1 /4 of the time Tgi- defined in 
equation ( |4TI) ; coalescence occurs in a time of Gyr if 

— « (0.015,0.034) X Ml ft^^^ (58) 
fl/, 

where the numbers in parentheses correspond to a = (1,0.1) 
respectively. The a/a/, values in equation ( |58] | are — 15 times 
smaller than those found by Sesana et al. ( 2007) after the sec- 
ondary slingshot had run its course, implying that the binaries 
in their model galaxies would stall at separations far outside 
th e gravitational radia tion regime unless extremely eccentric. 

Sesana et alj (l2007l) 's results might still be taken to imply 
that massive binaries commence their long-term, relaxation- 
driven evolution starting from separations somewhat smaller 
than ~ a/,, as assumed here. However such an effect was 
not apparent in the fully self-consistent A^-body simulations of 
iMerrittI (1200 6). This is probably due to the neglect by Sesana 
et al. of the changes in nuclear structure that accompany bi- 
nary formation. Sesana et al. computed the initial population 
of stars available to undergo reejections by assuming a singu- 
lar isothermal sphere density profile, p o= r^^, and counting 
the number of stars on orbits that intersected the binary. Even 
if such a steep density profile were present intially, it would 
be converted into a core of much lower density by the time 
a K, Gf,, and the number of stars available for the secondary 
slingshot would be much less than Sesana et al. estimated. 

In the Fokker-Planck integrations presented here (§5), the 
inclusion of the secondary slingshot had almost no effect on 
th e long-term behavior of a(f). 

IZied (12006 a.b ) also argued that stars near a binary black 
hole at the time of its formation could drive the binary to the 
gravitational radiation regime in a very short time. Zier ig- 
nored the secondary slingshot, but assumed that a dense clus- 
ter of stars would be bound to the binary at the time that its 
separation first reached ~ a/,. He found that a cluster having 
total mass ^ M12, distributed as a steep power-law around the 
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binary, p^r ^, Y<;2.5, could extract enough energy from it 
via the gravitational slingshot that Tgi- would fall below lO'" 
yr. While no detailed justification for such dense massive 
clusters was presented, Zier argued that "Low angular mo- 
mentum matter accumulates in the center" of merging galax- 
ies, and that "Each of the BHs will carry a stellar cusp with a 
mass of about its own" after the merger As noted above, re- 
cent observations do suggest the presence of compact nuclei 
at the centers of low-luminosity spheroids. 

A^-body simulatio ns of Zier's model hav e yet t o be car- 
ried out, although .Milosavljevic & Meirit^ (1200 I h did fol- 
low the evolution of merging galaxies with initial, p ~ 
r^^ cusps around each of the black holes, close to the 
value Y = 2.5 abo ve which Z ier infers rapid coalescence. 
iMilosavlievic & MerrittI (1200 lb observed a rapid phase of evo- 
lution of the binary, during which the density cusps were de- 
stroyed and a dropped by a factor of ~ a few below a/,. How- 
ever this was still far above the separation at which gravita- 
tional radiation would be efficient. 

An earl y, heuri stic model for binary evolution was pre- 
sented by MerrittI ( 12000) based on the results of the A^-body 
experiments that had been carried out up to that date. The 
model assumed that the rate of supply of stars to the binary 
was determined by local parameters (density, velocity disper- 
sion) and was independent of the nuclear relaxation time. The 
model was able to mimic the binary hardening rates seen in 
some A^-body experiments (Quinlan & Hernquist 1997; Chat- 
terjee, Hernquist & Loeb 2003); when scaled to real galax- 
ies, it predicted binary coalescence times that were nearly 
independent of galaxy mass. However, the A^-body results 
on which the model was based were subsequently called into 
question when they could not be reproduced using more ac- 
curate integrators (Makino & Funato 2004; Berczik, Merritt 
& S purzem 2005). M. Volonteri and co-authors adopted the 
iMerrittI (12000) prescription for binary evolution as a com- 
ponent of their semi-a n alytic models of black hole growth 
(IVolonterietani2003allR l2005h . and their models can there- 
fore be expected to substantially over-estimate the rate of bi- 
nary evolution in galaxies with M, > lO^M©. 

10. SUMMARY 

1. Accurate, long-term A^-body integrations of binary su- 
permassive black holes at the centers of realistically dense 
galaxy models were carried out using particle numbers up to 



0.26 X 10^. A new implementation of the Mikkola-Aarseth 
chain regularization algorithm was used to treat close interac- 
tions involving the black hole particles. The dependence of 
the binary's hardening rate on particle number was quantified 
by averaging the results of independent integrations. 

2. A Fokker-Planck model was developed that includes, for 
the first time, changes in the stellar density and potential due 
to star-binary interactions. The Fokker-Planck model was ver- 
ified by comparison with the averaged A^-body integrations. 

3. Based on the Fokker-Planck integrations and on empir- 
ical scaling relations, binary evolution in real galaxies was 
shown to take place in the "empty loss cone" (diffusive) 
regime for binaries with total mass above about This 
regime is out of range of particle numbers currently feasible 
via direct A^-body simulation but can be efficiently treated via 
the Fokker-Planck approximation. 

4. Accurate analytical expressions were derived that repro- 
duce the predictions of the Fokker-Planck model for the time- 
dependence of binary semi-major axis (equation [39]l and the 
time to coalescence (equation ISSTi in the diffusive regime. 

5. Based on the Fokker-Planck integrations and on em- 
pirical scaling relations, gravitational-radiation coalescence 
will occur in 10 Gyr or less for galaxies with binary masses 
< 2 X lO^'Mp; or central velocity dispersions 80 km s^^; 
the coalescence time depends only weakly on binary mass ra- 
tio (Fig. \T9[ . Binaries with masses ^ IO'^Mq will remam 
stalled for a Hubble time. 

6. A core, or "mass deficit," is created as a result of compe- 
tition between ejection of stars by the binary and re-supply of 
depleted orbits via gravitational (star-star) encounters. Mass 
deficits as large as ^ 4 (Mi +M2) were found to be generated 
before coalescence (Fig. ll6ll"7] i. 

7. After the black holes coalesce, a Bahcall-Wolf cusp 
forms around the single hole in approximately one relaxation 
time, resulting in a nuclear density profile with a flat core and 
an inner, compact cluster (Fig.l^TTi. similar to what is observed 
at the centers of low-luminosity spheroids. 
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